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Abstract — Using the spectral decomposition of the Laplace- 
Beltrami operator we simulate fractal surfaces as random series 
of eigenfunctions. This approach allows us to generate random 
fields over smooth manifolds of arbitrary dimension, generalizing 
previous work with fractional Brownian motion with multi- 
dimensional parameter. We give examples of surfaces with and 
without boundary and discuss implementation. 

Index Terms — Fractal surfaces, fractional Brownian motion, 
discrete Laplace-Beltrami operators. 

I. Introduction 

A Common approach to simulating fractal surfaces is via 
the sample paths of fractional Brownian motions and 
their multidimensional extensions to (e.g. O)- These 
random fields are self- similar in distribution in that when 
sampled at various scales the distribution of the sample is 
the same up to a constant scaling factor. However, they are 
indexed by . Suppose instead one wanted to simulate a 
fractal surface with some more complex geometry or topology, 
e.g., a fractal cylinder. The analogous approach to simulating 
such an object would require a random field indexed by a 
manifold, such as a surface in R^, that possessed the same 
properties as fractional Brownian motion, i.e., self- similarity, 
almost surely Holder continuous sample paths, and stationary 
increments. 

Constructing such fields has been a subject of ongoing 
research for a number of years and until recently, extensions 
were known only for limited classes of surfaces, such as 
spheres or hyperboloids, and only for certain ranges of self- 
similarity, i.e., only with Hurst index a G (0,1/2]. Then 
in |3 | self-similar Gaussian random fields were constructed 
over wide classes of surfaces, including arbitrary compact 
manifolds, both with and without boundary. These extensions 
were shown to have Holder continuous sample paths, be self- 
similar in distribution, and have stationary increments, thus 
possessing all the basic properties one would expect of an 
extension of fractional Brownian motion. In this paper we 
describe the simulation of these fields for surfaces in R^. 

The simulation of Gaussian random fields presents non- 
trivial challenges, in particular if one considers index sets 
other than R^. One common approach is the use of Fourier 
or wavelet techniques (e.g. ch. 2, |2|), but if one is 
working on a general surface the Fourier transform is no 
longer available. On the other hand, one can view the Fourier 
transform as the spectral decomposition of the Laplacian, the 
sine and cosine functions being the eigenfunctions. These we 
do have on a surface, and so the analogous approach is to 
build and simulate random fields over manifolds and surfaces 



using Fourier series of eigenfunctions of the Laplace-Beltrami 
operator. One novel feature of this approach is that using the 
Dirichlet Laplacian for a surface with boundary allows us to 
produce fields with almost sure boundary values. We begin 
the next section with a derivation in the simple case of the 
interval [0, 1], before moving on to surfaces. 

II. Fractional Brownian fields over surfaces 

A. A Basic Example 

As motivation and a base case for our construction, consider 
Brownian motion on [0,1]. Being a Gaussian process, Bt is 
determined by its covariance function, 

^[BtBs] min{s, t) = s M. 

Now notice that 5 A t is the Green's function of the Laplacian 
on [0, 1], -A - i.e., for / G L^p, 1] 

f^sMf{t)dt = f{s). (1) 

Of course we must specify boundary conditions to have a well 
defined Laplacian, and we see that these come from s f\t'. If 

F{s) = [ sMf{t)dt 
Jo 

then F(0) = and F\l) = 0. Thus we can say that Bt is 
the unique (up to equality in distribution) Gaussian process on 
[0, 1] whose covariance is the Green's function of —A acting 
on smooth functions / : [0, 1] ^ R such that /(O) = f{l) = 
0. 

If we let K : L^[0, 1] ^ L^[0, 1] be given by 

Kif){s)= I sMf{t)dt 
Jo 

then we can write (1) as 

K = {-A)-' (2) 

and in this way we see that starting from Bt, we uniquely 
determine — A as the inverse of the integral operator K defined 
by the covariance of Bt. If {A^}^^ and {0/c}^i are the 
eigenvalues and eigenfunctions respectively of —A, with the 
above boundary conditions, then a calculation shows 

Afc=^/c-0 TT^, 0/e(x) = V^sin ^/C - TTX. 

The Spectral Theorem and the functional calculus associated 
with it then yield 

CO 

-A(/) = ^Afe(/,,^fe)0fc, 
fe=i 
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Fig. 1. fractional Gaussian bridge with a = .1, .5, .9 from top to bottom. 
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and 



k=l 

(f^g) denoting the inner product, fgdx. If we now write 

K{f)= f\{x,y)f{y)dy, 
Jo 

oo 



k=l 



it is easily seen that the series defining k{x,y) converges 
absolutely and uniformly on [0, 1] and moreover it must be 
that k{x,y) = x Ay. Then if {^k} are i.i.d. standard normal 
random variables, an application of Fubini's Theorem yields 



E 



J2^khkMs)] (J2^khkMt) 



\k=l 



k{s,t) = sAt 



and we thus arrive at the well known Fourier series expansion 
of Bt, 



k=l 



k=l 



(3) 



Now let W be the white noise (or isonormal process, cf [?]) 
on L^[0, 1] and denote its action a function / G by 



W{f) ^ J 



fdW. 



Then {W{(})k)} is a set of i.i.d. standard normal random 
variables, denoted again by {ik}^ and if 

oo 

k^{x,y) = 2_^ \ ^ (l)k{x)(l)k{y) 

.{{k- ^)7Tx) sin {{k - I) Try) 



k=l 



(fc-l)7r 



then 



{-A)-- f{x)= / kHx,y)f{y)dy 
Jo 



and 



(-A)-M^) = I kHx,y)dW{y) 

oo ^ 



k=l 



sin (/c — ^) irt 



_ 1 

Thus we can write Bt — (—A) {W) and arrive at the 
stochastic steady- state equation 



{-AYBt = W. 



(4) 



Notice that working backwards we can define Bt to be the 
unique Gaussian process satisfying (4). This is similar, but 
not the same as, the classical equation 



~dt 



Bt = W. 



To see that (4) is in fact different, note that (— A)"^ is self 
adjoint on L^[0, 1], while d/dt is not (integrate by parts). 

Now suppose we replace the mixed Laplacian above by the 
Dirichlet Laplacian in (4), acting on functions such that /(O) = 
/(I) = 0. The Dirichlet Laplacian on [0, 1] has eigenvalues 
and eigenfunctions given by 



A. 



(kir)^ 



and 



V2sm{k7rx) 

respectively. Then we have as above a process Xt defined by 

sin(/c7rt) 



k=i 



kiT 



If we now let a G (0, 1) (the case above being a = 1/2) we 
can extend (4) to 



and write 



(_A)(^ + t)x^ = W 



— 2_^^k /I , V 



(5) 



(6) 



which is the homogeneous Riesz field constructed in |3|. This 
is a self-similar Gaussian process with almost surely Holder 
continuous sample paths and such that Xq = Xi = almost 
surely, i.e., it is a Gaussian bridge. 

Because the series converges in with probability 1, we 
can simulate these processes by taking partial sums (see fig. 
1). Notice that as a increases the sample paths become more 
regular. The reason can be seen from (6): Lower alpha em- 
phasizes the larger eigenvalues and thus the higher frequency 
wave functions show more of a contribution. Similarly a larger 
a suppresses their contribution and the lower frequency wave 
functions dominate, leading to less overall oscillation and a 
smoother sample path. We will see this same behavior in 
general below. 
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(C) a = .9 

Fig. 2. Cylinders with Dirichlet boundary conditions for o; = .1, .5, .9 

B. Extension to Surfaces with Boundary 

Suppose now D is a. compact connected surface in 
with smooth boundary. Then there is a canonical differential 
operator on C^{D) corresponding to the Euclidean Laplace 
operator above, the Laplace-Beltrami operator, which we refer 
to also as simple the Laplacian of D (see e.g. |4|). As above 
we let \k and (j)k be the eigenvalues and eigenfunctions of the 
Dirichlet Laplacian of D, —A. Given a white noise W on D 
we start from 

(-A)(^ + ^)i?^ = W 

to obtain 

/c=l 



where the change in the exponent from 1 to 2 ensures sample 
path regularity (see |3 1). In this way we obtain a random field 
on D that is self-similar, almost-surely Holder continuous of 
order 7 for any 7 < a and almost surely not Holder continuous 
for 7 > a, and has stationary increments. A few words are 
in order as to what we mean by self-similar: Because we are 
considering surfaces embedded in R^, that is, surfaces with 
Riemannian metrics induced by the Euclidean metric on R^, 
the usual definitions carry over, i.e.. Rex = c^Rx- For more 
general considerations of self- similarity on arbitrary manifolds 
see 0. 

Now we have a generalization of the bridges above that 
possesses the properties we would like, however for a general 
surface D we no longer have simple analytical expressions 
for the functions (j)k. One way around this is to discretize A 
and construct the analogous discretized version of (7). This 
requires some care, as there are a number of discretizations 
available for Laplace-Beltrami operators, some of which be- 
have very differently (e.g. |5|). 

Given a meshed surface, most discrete Laplacians are 
weighted graph Laplacians, 

Lf{x) = Y,w,,y{f{x)-f{y)) 

for some edge weights Wx^y and where x y means x and y 
are adjacent in the mesh (i.e. they share an edge). Supposing 
we have a fairly uniform mesh, let us SQi Wx^y = ||x — ?/||^3^. 
Denote this discrete Laplacian by L, and let its eigenvectors 
and eigenvalues be {(^/e} and {\k} respectively. Note that ^ 
R^ where the dimension of L (L is a matrix) is N x N, and 
is a real valued function on our mesh. We then let 

Rx = Y.~^k^ iM^). (8) 

k=l 

where we choose some Nq < N and ^/e are again independent 
standard normals. 

We would like to be able to take our mesh fine enough and 
A^o large enough so that Rx provides a good approximation 
of the full Riesz field Rx. How good this approximation is, 
and in what sense, will depend completely on how well L 
approximates A. For general surfaces, discrete Laplacians and 
their convergence to the continuum Laplacian is a very active 
area of current research, and we defer a more full discussion 
for the moment. 

Some examples of planar domains are plotted in Figs. 3, 
5, and 6. Note that with probability 1 the field is zero on 
the boundary. As an example of what is possible in three 
dimensions, consider the cylinder [0, C] x (Fig. 2). If 
we impose zero Dirichlet boundary conditions on the two 
boundary circles, we get analogues of the bridges in Fig. L 
Notice that again, a controls the sample path (or surface) 
regularity and the fields are almost surely equal to zero on 
both boundary circles, i.e., the height of the field values plotted 
along normals is zero. 

C. Surfaces Without Boundary 

For surfaces without boundary there are no obvious bound- 
ary conditions, so we need to interpret our equations differ- 



(a) a = .1 




Fig. 3. Disk with 3 smaller disks removed, a = .9 

ently. The basic issue is that on a compact manifold, e.g., 
the sphere S^, the lowest eigenvalue of the Laplacian is zero 
and (— A)~^2 + f ) is not defined. We can instead consider the 
series 



/c=l 



(9) 



where o is some fixed point serving as an "origin." This 
series does converge and again defines a self- similar field with 
stationary increments and almost sure Holder continuity as 
above, called the Riesz Field in |3|. Figure 4 shows spheres 
with increasing alpha, using the same discrete Laplacian L as 
before. 




(b) a = .5 




Fig. 4. Spheres with a ■ 
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III. Implementation 

A. Convergence 

For uniform or nearly uniform meshes, the above discretiza- 
tion produces good results. In general however, placing a 
uniform mesh on a surface is non-trivial (depending on the 
curvature and boundary), and in this case the above discrete 
Laplacian may be a poor approximation to the continuous 
Laplacian. In the case of highly curved surfaces, care must 
be taken in choosing an appropriate discretization. 

Suppose we have a meshed surface and have chosen a 
discrete Laplacian, L. A simple condition for the almost sure 
approximation — ^||x,2 < e is as follows: Suppose that 
for any integer A/" > 0, for all n < we have 

lim - 0n||L2 

mesh^o 



and 



|An - \n\ 0, 



where we can embed in L'^{D^dV), dV denoting surface 
measure on D. Then for any e > we can choose A'e such 
that with probability 1 



lim 
mesh^ 



Ne 



n=l 

< e. 



oo 
n=l 



-(f + f) 



L2 



This is a simple consequence of the fact that 
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converges in LP' almost surely, for then we just bound the tail 
and the assumed convergence of the lower eigenfunctions and 
eigenvalues yields the claim. 

Thus if one has a meshed surface, the Laplacian chosen 
should satisfy the above spectral convergence condition. Such 
conditions are known to hold for certain Laplacians, but not 
all I61, 13 . Of course, this is an active area of research and 
so we expect further results on spectral approximation of the 
Laplace-Beltrami operator to become applicable to the spectral 
synthesis we have used here. 

B. Efficiency and Self -Similarity 

There are two parameters which determine the computation 
time in simulating for a given surface: The number n of 
grid points in the mesh (i.e. how fine the mesh is) and the 
number N of eigenvectors appearing in 

N (1 \ 

1 

The dimension of the discrete Laplacian will be n x n and in 
general N should be taken much smaller than n^, for Weyl's 
asymptotic formula for the eigenvalues of A (e.g. |4|) tells 
us that Xk = 0{k) for a two dimensional surface as we are 
dealing with here. Thus 

A-(^+^)=o(fc-(i+f)) 

and so if our discrete Laplacian L satisfies the spectral 
convergence above, we see that the contribution, in mean 
square, of the A^th term is of the order A^~(2+2^)^ which 
may rapidly become negligible depending on the fineness of 
the mesh. 

On the other hand, for finer meshes, the computation of 
even the first few hundred eigenvectors and eigenvalues of the 
n X n matrix L can become long. However, once computed, 
this spectral data can be stored for repeated use. All that is 
needed is the random numbers {^/c}, which are not expensive. 
Moreover, the self- similarity of the Riesz fields allows one to 
compute the spectra of a surface of any size, and simply rescale 
the resulting field without computing the spectra for a rescaled 
surface. For example, once one has simulated the Riesz field 
on a sphere of radius 1, one can simulate it on a sphere of 
radius c simply by plotting c^Rx for the appropriate value of 
a. In this way, for very detailed images, the computational 
time required may be somewhat high, but it is a one time cost 
in the above sense. 

IV. Possible Extensions 

We have restricted our attention to embedded surfaces in R^, 
however the theory in 1 3 1 applies to non-embedded manifolds, 
e.g., hyperbolic space or flat tori, as well as n-dimensional 
surfaces in R^+^. In this case the approximation takes the 
form 

1 

Of course working with a surface whose metric is not induced 
by the ambient Euclidean metric increases the difficulty in 



Fig. 5. Fractal surface on the unit disk with Dirichlet boundary conditions 
{a = .7) with negative values set to zero. 




(b) 

Fig. 6. Fractal surface over a teardrop shaped region with a = .7. 



computing accurate discretizations of the Laplacian. Also, we 
have here chosen a to be fixed, but one could also allow a to 
vary over the surface and thus obtain multifractal fields. We 
expect, similar to what is known on R^, that such fields with 
varying fractal index would still be continuous. 
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